Multivariate TS Models (ARIMAX/SARIMAX/VAR)

Literature Review

In analyzing the temporal dynamics of the U.S. freight transportation industry, key variables emerge as critical factors shaping its trajectory. Economic indicators, encompassing metrics like GDP growth, consumer spending, and industrial production, serve as pivotal markers of overall economic activity. These indicators not only influence freight transportation demand but also offer insights into broader economic trends. Additionally, fuel prices wield significant influence, affecting operating costs and modal choice decisions within the transportation sector. Research underscores the impact of fuel price fluctuations on transportation costs, suggesting their inclusion as exogenous variables in ARIMAX models to enhance freight volume forecasts.

When considering modeling approaches, both ARIMAX and VAR models offer distinct advantages for capturing the complexities of freight transportation dynamics. ARIMAX models are well-suited for incorporating economic indicators, fuel prices, and infrastructure investment as variables. Leveraging their autoregressive and exogenous components, ARIMAX models can effectively account for temporal patterns and external influences, thereby improving the accuracy of freight volume forecasts over time. On the other hand, VAR models excel in analyzing the dynamic relationships among multiple variables. They provide a versatile framework for examining the simultaneous effects of economic indicators, fuel prices, regulatory changes, and technological advancements on freight transportation dynamics. With their capacity to capture intricate interactions and feedback mechanisms, VAR models offer valuable insights into the evolving landscape of the freight transportation industry.

In my ARIMAX/SARIMAX model, I’ve chosen to focus on two key datasets. Firstly, the COVID-19 pandemic has reshaped societal norms significantly, leading to notable shifts in consumer behaviors and lifestyle patterns. With restrictions in place and mobility constrained, there has been a substantial surge in demand for delivery services. UPS, as a prominent player in the U.S. delivery service sector, has been instrumental in meeting this increased demand. My goal is to examine whether there exists a relationship between the number of new COVID-19 cases and UPS stock prices. To accomplish this, I will employ an ARIMAX model to delve into the dynamics between these variables and discern any potential impact of the pandemic on UPS stock performance.

Secondly, fuel prices play a pivotal role in the freight industry, exerting a significant influence on operational costs. I aim to investigate whether fluctuations in fuel prices affect freight values, particularly within the international freight sector. Additionally, I aim to assess the potential impact of the Transportation Services Index (TSI) on freight values. To explore these relationships comprehensively, I intend to utilize the SARIMAX model. This model facilitates the analysis of time series data while accounting for external factors such as fuel prices and the TSI. Through this analysis, I seek to gain insights into the interplay between fuel prices, the TSI, and the values of Canadian freight in the U.S. landscape.

For the VAR model, the analysis delves into three distinct sets of variables within the air cargo industry to uncover their interrelationships and implications. Firstly, it examines the connections between Air Cargo (Domestic), Air Cargo (International), Fuel Price, and the Transportation Services Index (TSI). This investigation aims to discern whether fluctuations in air cargo volumes, both domestically and internationally, correlate with changes in fuel prices and TSI, thereby providing insights into the broader drivers of the air cargo industry.

Secondly, the analysis focuses on the relationship between Air Cargo (Domestic), Air Cargo (International), and Employment. By exploring the mutual influence between air cargo volumes and employment levels. Understanding these dynamics is crucial for grasping the industry’s labor market dynamics and its broader economic implications.

Lastly, the analysis considers the dynamics between GDP, Revenue, and Employment within the air industry. By examining the interplay between these key economic indicators, it aims to uncover the industry’s economic performance and its effects on employment trends. This exploration provides valuable insights into the industry’s contribution to the broader economy and its role in shaping employment patterns.

Define Models

SARIMAX: UPS Stock Price ~ Covid New Cases (Weekly)

Code
#  SARIMAX: ups stock price and covid new cases (weekly)
# get the weekly stock price of ups
ups <- getSymbols("UPS",auto.assign = FALSE, from = "2020-01-01", to = "2023-05-15",src="yahoo") 
ups=data.frame(ups)
ups <- data.frame(ups,rownames(ups))
colnames(ups)[7] = "date"
ups$date<-as.Date(ups$date,"%Y-%m-%d")
ups <- ups[seq(1, nrow(ups), by = 5), ]

# get the weekly covid new cases
covid <- read.csv("data/clean_data/covid.csv",header =TRUE, sep = ',')
covid <- subset(covid, select = c(date, new_cases_per_million))
covid <- na.omit(covid)
covid <- head(covid, -6)
covid$new_cases_per_million[covid$new_cases_per_million == 0.000] <- 0.001

# combine the data
ups_covid <- data.frame("UPS"=log(ups$UPS.Adjusted), "COVID"=log(covid$new_cases_per_million))
#(head(ups_covid))

# Convert to a Time Series Component 
ups_covid_ts <- ts(ups_covid, start = decimal_date(as.Date("2020-01-02", format = "%Y-%m-%d")), frequency = 52)

head(ups_covid_ts)
Time Series:
Start = 2020.00273224044 
End = 2020.09888608659 
Frequency = 52 
              UPS     COVID
2020.003 4.623659 -6.907755
2020.022 4.624857 -6.907755
2020.041 4.639122 -6.907755
2020.060 4.622545 -3.863233
2020.080 4.503046 -5.115996
2020.099 4.499175 -3.411248
Code
#VARselect(ups_covid_ts, lag.max=10, type="both")
#summary(vars::VAR(ups_covid_ts, p=3, type='both'))

SARIMAX: Canada Freight Value ~ Fuel Price + TSI (Monthly)

Code
#  SARIMAX: Canada Freight Value and Fuel Price, TSI (monthly)
# get the canada freight monthly value
canada <- read.csv("data/clean_data/24monthly_canada_freight.csv")
canada <- canada[,c("Date","Value")]

# get the fuel monthly price
fuel <- read.csv("data/clean_data/34monthly_fuel_prices.csv")
# remove the dollar sign
fuel$Diesel <- as.numeric(gsub("\\$", "", fuel$Diesel))
fuel$Jet.Fuel <- as.numeric(gsub("\\$", "", fuel$Jet.Fuel))
# select the same start data
fuel1 <- fuel[73:(nrow(fuel)-6),]

# get the TSI monthly data
tsi <- read.csv("data/clean_data/35monthly_TSI.csv")
tsi1 <- tsi[73:nrow(tsi),]

# combine the data
canada_value <- data.frame("Canada_value"=canada$Value, "Fuel"=fuel1$Diesel, "TSI"=tsi1$Transportation.Services.Index...Freight)
#(head(canada_value))

canada_value_ts <- ts(canada_value, start = decimal_date(as.Date("2006-01-01", format = "%Y-%m-%d")), frequency = 12)

head(canada_value_ts)
         Canada_value  Fuel   TSI
Jan 2006     43.47682 246.7 112.7
Feb 2006     42.49097 247.5 111.8
Mar 2006     47.70424 255.9 112.3
Apr 2006     44.39409 272.8 111.5
May 2006     46.93539 289.7 113.8
Jun 2006     46.93349 289.8 113.2
Code
#VARselect(canada_value_ts, lag.max=10, type="both")
#summary(vars::VAR(canada_value_ts, p=2, type='both'))

VAR: Air Cargo(Domestic) ~ Air Cargo(International) + Fuel Price + TSI (Monthly)

Code
#  VAR: Air Cargo and Fuel Price, TSI (monthly)
# get the cargo data
cargo <- read.csv("data/clean_data/Monthly_air_revenue.csv")
cargo$International <- cargo$U.S..Air.Carrier.Cargo..millions.of.revenue.ton.miles....International/1000000000
cargo$Domestic <- cargo$U.S..Air.Carrier.Cargo..millions.of.revenue.ton.miles....Domestic/1000000000
cargo <- cargo[4:(nrow(cargo)-8),]

# get the fuel monthly price
fuel <- read.csv("data/clean_data/34monthly_fuel_prices.csv")
# remove the dollar sign
fuel$Diesel <- as.numeric(gsub("\\$", "", fuel$Diesel))
fuel$Jet.Fuel <- as.numeric(gsub("\\$", "", fuel$Jet.Fuel))
# select the same start data
fuel2 <- fuel[37:(nrow(fuel)-6),]

# get the TSI monthly data
tsi <- read.csv("data/clean_data/35monthly_TSI.csv")
tsi2 <- tsi[37:nrow(tsi),]

# combine the data
revenue_ <- data.frame("Domestic"=cargo$Domestic, "International"=cargo$International, "Fuel"=fuel2$Jet.Fuel, "TSI"=tsi2$Transportation.Services.Index...Freight)
#(head(revenue_))

revenue_ts <- ts(revenue_, start = decimal_date(as.Date("2003-01-01", format = "%Y-%m-%d")), frequency = 12)

head(revenue_ts)
         Domestic International  Fuel   TSI
Jan 2003 1.177785      1.395876  88.7 104.2
Feb 2003 1.079561      1.366617 105.5 103.4
Mar 2003 1.233953      1.691870  89.3 103.9
Apr 2003 1.208202      1.524837  74.3 104.0
May 2003 1.243021      1.489449  71.4 103.0
Jun 2003 1.207785      1.473214  74.8 102.9

VAR: Employment ~ Air Cargo(Domestic) + Air Cargo(International) (Monthly)

Code
#  VAR: Employment and Air Cargo (monthly)
# get the employment data
employment <- read.csv("data/clean_data/32monthly_employment.csv")
employment <- employment[,c("Date","Transportation.Employment...Air.Transportation")]

# get the cargo data
cargo <- read.csv("data/clean_data/Monthly_air_revenue.csv")
cargo$International <- cargo$U.S..Air.Carrier.Cargo..millions.of.revenue.ton.miles....International/1000000000
cargo$Domestic <- cargo$U.S..Air.Carrier.Cargo..millions.of.revenue.ton.miles....Domestic/1000000000
cargo1 <- cargo[28:(nrow(cargo)-8),]

# combine the data
emp_ <- data.frame("Employment"=employment$Transportation.Employment...Air.Transportation, "International"=cargo1$International, "Domestic"=cargo1$Domestic)
#(head(emp_))

emp_ts <- ts(emp_, start = decimal_date(as.Date("2005-01-01", format = "%Y-%m-%d")), frequency = 12)

head(emp_ts)
         Employment International Domestic
Jan 2005     505200      1.696944 1.280516
Feb 2005     503400      1.561959 1.253130
Mar 2005     504200      1.958162 1.469765
Apr 2005     507300      2.078222 1.226824
May 2005     507800      1.925489 1.261165
Jun 2005     508300      1.981247 1.351714

VAR: GDP ~ Revenue + Employment (Air, Yearly)

Code
#  VAR: Air GDP, Revenue and Employment (yearly)
# get the yearly GDP of air
gdp <- read.csv("data/clean_data/31GDP.csv")
# Group by Year and Mode, calculate the sum of GDP within each group
gdp <- gdp %>%
  filter(Mode =='Air')%>%
  group_by(Year) %>%
  summarize(Value = sum(Current....billions., na.rm = TRUE))

# get the air revenue
air <- read.csv("data/clean_data/33revenue.csv")
air <- air[air$Mode=="Air carrier, domestic",c("Year","Value")]
air <- air[order(air$Year), ]
air <- air[13:nrow(air),]

# get the air employment
employment <- read.csv("data/clean_data/32monthly_employment.csv")
employment <- employment[,c("Date","Transportation.Employment...Air.Transportation")]
employment$Date <- as.POSIXct(employment$Date, format = "%m/%d/%Y %I:%M:%S %p")
employment$Year <- format(employment$Date, "%Y")
employment <- employment %>%
              group_by(Year) %>%
              summarize(Value = sum(Transportation.Employment...Air.Transportation, na.rm = TRUE)/1000000)
employment <- employment[8:17,]

# combine the data
air_ <- data.frame("GDP"=gdp$Value, "Revenue"=air$Value, "Employment"=employment$Value)
#(head(air_))

air_ts <- ts(air_, start = decimal_date(as.Date("2012-01-01", format = "%Y-%m-%d")), frequency = 1)

(air_ts)
Time Series:
Start = 2012 
End = 2021 
Frequency = 1 
           GDP  Revenue Employment
2012 115.38111 133.5094     5.5101
2013 122.50650 135.4657     5.3318
2014 127.33886 135.8515     5.3308
2015 139.52287 124.0370     5.5027
2016 152.46236 119.8993     5.7305
2017 164.88538 125.1505     5.8978
2018 174.17990 137.4876     5.9631
2019 183.64735 137.6700     6.0585
2020  74.96704 122.9800     5.2947
2021 105.60998  97.8600     5.5136

ARIMAX/SARIMAX

Select Variables

Code
#ups_covid_ts
autoplot(ups_covid_ts, facets = T) +
  labs(x="Date", y="", title = "UPS Stock Price and Covid New Cases")

The COVID-19 pandemic has been a profound influence on society, dramatically altering people’s lifestyles and behaviors. With restrictions in place and mobility limited, there has been a significant surge in demand for delivery services. UPS, as a leading company in the delivery service sector in the U.S., stands at the forefront of meeting this increased demand. I aim to investigate whether there is a relationship between the number of new COVID-19 cases and UPS stock prices. To do this, I will utilize an ARIMAX model to explore the dynamics between these variables and discern any potential impact of the pandemic on UPS stock performance.

Code
#canada_value_ts
autoplot(canada_value_ts, facets = T) +
  labs(x="Date", y="", title = "Canada Freight Value and Fuel Price and TSI")

Fuel prices play a crucial role in the freight industry, significantly impacting operational costs. I aim to investigate whether fluctuations in fuel prices influence freight values, particularly in the international freight sector. Additionally, I seek to understand the potential impact of the Transportation Services Index (TSI) on freight values. To explore these relationships, I plan to utilize the SARIMAX model, which allows for the analysis of time series data while considering exogenous variables such as fuel prices and the TSI. By examining these dynamics, I hope to gain insights into the interplay between fuel prices, the TSI, and the Canada freight values in the U.S..

Auto.Arima()

Code
#ups_covid_ts
fit0 <- auto.arima(ups_covid_ts[,'UPS'], xreg=ups_covid_ts[,'COVID'])
summary(fit0)
Series: ups_covid_ts[, "UPS"] 
Regression with ARIMA(0,1,0)(1,0,0)[52] errors 

Coefficients:
         sar1     xreg
      -0.0795  -0.0008
s.e.   0.0930   0.0068

sigma^2 = 0.001926:  log likelihood = 289.38
AIC=-572.75   AICc=-572.61   BIC=-563.36

Training set error measures:
                      ME       RMSE        MAE        MPE      MAPE      MASE
Training set 0.003119728 0.04349475 0.03097627 0.06025101 0.6193865 0.1143339
                   ACF1
Training set 0.01248485
Code
checkresiduals(fit0)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(0,1,0)(1,0,0)[52] errors
Q* = 32.331, df = 33, p-value = 0.5002

Model df: 1.   Total lags used: 34

Select model: Regression with ARIMA(0,1,0)(1,0,0)[52] errors.

Code
#canada_value_ts
xreg <- cbind(fuel = canada_value_ts[, "Fuel"],
              tsi = canada_value_ts[, "TSI"])
fit0 <- auto.arima(canada_value_ts[,'Canada_value'], xreg=xreg)
summary(fit0)
Series: canada_value_ts[, "Canada_value"] 
Regression with ARIMA(1,0,0)(1,0,0)[12] errors 

Coefficients:
         ar1    sar1  intercept    fuel     tsi
      0.6281  0.6640   -23.8772  0.0581  0.4504
s.e.  0.0581  0.0602     8.0465  0.0059  0.0726

sigma^2 = 5.295:  log likelihood = -462.94
AIC=937.88   AICc=938.3   BIC=957.81

Training set error measures:
                      ME     RMSE      MAE        MPE     MAPE    MASE
Training set -0.05212728 2.272867 1.762864 -0.2909279 3.694124 0.30965
                   ACF1
Training set -0.0597277
Code
checkresiduals(fit0)


    Ljung-Box test

data:  Residuals from Regression with ARIMA(1,0,0)(1,0,0)[12] errors
Q* = 42.707, df = 22, p-value = 0.005127

Model df: 2.   Total lags used: 24

Select model: Regression with ARIMA(1,0,0)(1,0,0)[12] errors.

Manual Method

Code
#ups_covid_ts
fit.reg1 <- lm(UPS ~ COVID, data=ups_covid_ts)
summary(fit.reg1)

Call:
lm(formula = UPS ~ COVID, data = ups_covid_ts)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.62194 -0.08510  0.06658  0.13999  0.29465 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.70663    0.04087 115.157  < 2e-16 ***
COVID        0.04988    0.00577   8.645 4.11e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2164 on 168 degrees of freedom
Multiple R-squared:  0.3079,    Adjusted R-squared:  0.3038 
F-statistic: 74.73 on 1 and 168 DF,  p-value: 4.114e-15

COVID is significant.

Code
res.fit1<-ts(residuals(fit.reg1),star=decimal_date(as.Date("2020-01-02",format = "%Y-%m-%d")),frequency = 52)

ggtsdisplay(res.fit1)

The data is not stationary, it needs differencing.

Code
ggtsdisplay(diff(res.fit1))

p=2,3,d=1,q=2,3,P=1,D=0,Q=0

Code
SARIMA.c=function(p1,p2,q1,q2,P1,P2,Q1,Q2,data){
  
  temp=c()
  d=0
  D=0
  s=52
  
  i=1
  temp= data.frame()
  ls=matrix(rep(NA,9*150),nrow=150)
  
  
  for (p in p1:p2)
  {
    for(q in q1:q2)
    {
      for(P in P1:P2)
      {
        for(Q in Q1:Q2)
        {
          if(p+d+q+P+D+Q<=10)
          {
            
            model<- Arima(data,order=c(p-1,d,q-1),seasonal=c(P-1,D,Q-1))
            ls[i,]= c(p-1,d,q-1,P-1,D,Q-1,model$aic,model$bic,model$aicc)
            i=i+1
            #print(i)
            
          }
          
        }
      }
    }
    
  }
  
  temp= as.data.frame(ls)
  names(temp)= c("p","d","q","P","D","Q","AIC","BIC","AICc")
  
  temp
  
}

# p=2,3,d=1,q=2,3,P=1,D=0,Q=0

output=SARIMA.c(p1=1,p2=4,q1=1,q2=3,P1=1,P2=1,Q1=1,Q2=1,data=res.fit1)
#output

output[which.min(output$AIC),] 
   p d q P D Q       AIC       BIC      AICc
12 3 0 2 0 0 0 -535.9344 -513.9838 -535.2431
Code
output[which.min(output$BIC),]
  p d q P D Q       AIC       BIC      AICc
8 2 0 1 0 0 0 -534.3201 -518.6411 -533.9542
Code
output[which.min(output$AICc),]
   p d q P D Q       AIC       BIC      AICc
12 3 0 2 0 0 0 -535.9344 -513.9838 -535.2431

Best models: ARIMA(3,0,2)x(0,0,0)[52], auto.arima suggested: ARIMA(0,1,0)(1,0,0)[52]

Model diagnostics to compare the two models.

Model1: ARIMA(3,0,2)x(0,0,0)[52]

Code
# model1: ARIMA(3,0,2)x(0,0,0)[52]
set.seed(1234)
model_output123 <- capture.output(sarima(res.fit1, 3,0,2, 0,0,0,52)) 

Code
cat(model_output123[105:125], model_output123[length(model_output123)], sep = "\n")
iter  44 value -3.036392
iter  45 value -3.036392
iter  46 value -3.036392
iter  47 value -3.036393
iter  48 value -3.036393
iter  48 value -3.036393
iter  48 value -3.036393
final  value -3.036393 
converged
<><><><><><><><><><><><><><>
 
Coefficients: 
      Estimate     SE t.value p.value
ar1     1.0032 0.1123  8.9331  0.0000
ar2     0.7658 0.1589  4.8208  0.0000
ar3    -0.7919 0.0985 -8.0412  0.0000
ma1     0.0095 0.1403  0.0675  0.9463
ma2    -0.6163 0.1250 -4.9310  0.0000
xmean   0.0177 0.0618  0.2867  0.7747

sigma^2 estimated as 0.002259191 on 164 degrees of freedom 
 

Model2: ARIMA(0,1,0)(1,0,0)[52]

Code
# model2: ARIMA(0,1,0)(1,0,0)[52]
set.seed(1234)
model_output12 <- capture.output(sarima(res.fit1, 0,1,0, 1,0,0,52)) 

Code
cat(model_output12[20:27], model_output12[length(model_output12)], sep = "\n")
Coefficients: 
         Estimate     SE t.value p.value
sar1      -0.1573 0.0945 -1.6643  0.0979
constant  -0.0008 0.0035 -0.2249  0.8223

sigma^2 estimated as 0.00250965 on 167 degrees of freedom 
 
AIC = -3.106522  AICc = -3.106094  BIC = -3.050962 
 

The Standard Residual Plot both appears to have some irregularities in the mean and variation.

The Autocorrelation Function (ACF) plot model1 shows almost no correlation indicating that the model has harnessed everything and all that is left is white noise. This indicates a good model fit. model2 shows some correlations left.

The Quantile-Quantile (Q-Q) Plot both shows near normality.

The Ljung-Box test results model1 reveal values above the 0.05 (5% significance) threshold, indicating a good fit.

$ttable: All coefficients are significant below a 0.05 threshold except one in model1.

When comparing the two model diagnostics, the manual model is better, which was Regression with ARIMA(3,0,2)(0,0,0)[52] errors.

Code
#canada_value_ts
fit.reg <- lm( Canada_value ~ Fuel+ TSI, data=canada_value_ts)
summary(fit.reg)

Call:
lm(formula = Canada_value ~ Fuel + TSI, data = canada_value_ts)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.3748  -2.0957   0.2293   2.0165   6.5057 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3.963655   2.242470  -1.768   0.0786 .  
Fuel         0.067373   0.002965  22.726   <2e-16 ***
TSI          0.260343   0.018099  14.384   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.143 on 202 degrees of freedom
Multiple R-squared:  0.8101,    Adjusted R-squared:  0.8082 
F-statistic: 430.8 on 2 and 202 DF,  p-value: < 2.2e-16

Fuel and TSI both be significant.

Code
res.fit<-ts(residuals(fit.reg),star=decimal_date(as.Date("2006-01-01",format = "%Y-%m-%d")),frequency = 12)

ggtsdisplay(res.fit)

The data is not stationary, it needs differencing.

Code
ggtsdisplay(diff(res.fit, lag=12))

p=1,d=0,q=1,2,3,P=1,2,D=1,Q=1

Code
SARIMA.c=function(p1,p2,q1,q2,P1,P2,Q1,Q2,d1,d2,data){
  
  temp=c()
  d=0
  D=1
  s=12
  
  i=1
  temp= data.frame()
  ls=matrix(rep(NA,9*150),nrow=150)
  
  
  for (p in p1:p2)
  {
    for(q in q1:q2)
    {
      for(P in P1:P2)
      {
        for(Q in Q1:Q2)
        {
          for(d in d1:d2)
       
        {
          if(p+d+q+P+D+Q<=10)
          {
            
            model<- Arima(data,order=c(p-1,d,q-1),seasonal=c(P-1,D,Q-1))
            ls[i,]= c(p-1,d,q-1,P-1,D,Q-1,model$aic,model$bic,model$aicc)
            i=i+1
            #print(i)
            
          }
          
        }
      }
    }
    
  }
  
  }
  temp= as.data.frame(ls)
  names(temp)= c("p","d","q","P","D","Q","AIC","BIC","AICc")
  
  temp
  
}

# p=1,d=0,q=1,2,3,P=1,2,D=1,Q=1

output=SARIMA.c(p1=1,p2=2,q1=1,q2=2,P1=1,P2=3,Q1=1,Q2=2,d1=0,d2=1,data=res.fit)
#output

output[which.min(output$AIC),] 
   p d q P D Q    AIC      BIC     AICc
35 1 0 0 2 1 1 849.28 865.5935 849.6009
Code
output[which.min(output$BIC),]
   p d q P D Q    AIC      BIC     AICc
35 1 0 0 2 1 1 849.28 865.5935 849.6009
Code
output[which.min(output$AICc),]
   p d q P D Q    AIC      BIC     AICc
35 1 0 0 2 1 1 849.28 865.5935 849.6009

Best models: ARIMA(1,0,0)x(2,1,1)[12], auto.arima suggested: ARIMA(1,0,0)(1,0,0)[12]

Model diagnostics to compare the two models.

Model1: ARIMA(1,0,0)x(2,1,1)[12]

Code
# model1: ARIMA(1,0,0)x(2,1,1)[12]
set.seed(1234)
model_output12 <- capture.output(sarima(res.fit, 1,0,0, 2,1,1,12)) 

Code
cat(model_output12[30:40], model_output12[length(model_output12)], sep = "\n")
Coefficients: 
         Estimate     SE t.value p.value
ar1        0.6695 0.0542 12.3449  0.0000
sar1       0.1933 0.0873  2.2151  0.0280
sar2      -0.2257 0.0843 -2.6785  0.0081
sma1      -0.9132 0.1159 -7.8776  0.0000
constant  -0.0059 0.0075 -0.7868  0.4324

sigma^2 estimated as 4.001151 on 188 degrees of freedom 
 
AIC = 4.407641  AICc = 4.409304  BIC = 4.509072 
 

Model2: ARIMA(1,0,0)x(1,0,0)[12]

Code
# model2: ARIMA(1,0,0)x(1,0,0)[12]
set.seed(1234)
model_output12 <- capture.output(sarima(res.fit, 1,0,0, 1,0,0,12)) 

Code
cat(model_output12[41:49], model_output12[length(model_output12)], sep = "\n")
Coefficients: 
      Estimate     SE t.value p.value
ar1     0.5987 0.0575 10.4168  0.0000
sar1    0.5922 0.0594  9.9712  0.0000
xmean   0.0671 0.9137  0.0734  0.9415

sigma^2 estimated as 5.427344 on 202 degrees of freedom 
 
AIC = 4.595809  AICc = 4.596392  BIC = 4.660648 
 

The Standard Residual Plot both appears to have some irregularities in the mean and variation.

The Autocorrelation Function (ACF) plot model1 shows almost no correlation indicating that the model has harnessed everything and all that is left is white noise. This indicates a good model fit. model2 shows some correlations left.

The Quantile-Quantile (Q-Q) Plot both shows near normality.

The Ljung-Box test results model1 reveal values above the 0.05 (5% significance) threshold, indicating a good fit.

$ttable: All coefficients are significant below a 0.05 threshold.

When comparing the two model diagnostics, the manual model is better, which was Regression with ARIMA(1,0,0)(2,1,1)[12] errors.

Cross-Validation

Code
#ups_covid_ts
n=length(res.fit1) #170  156
k=52

 #n-k=104; 104/52=2;

res <- res.fit1[1:156] 
rmse1 <- matrix(NA, 2,52)
rmse2 <- matrix(NA, 2,52)

st <- tsp(res.fit1)[1]+(k-1)/52 

for(i in 1:2)
{
  xtrain <- window(res.fit1, end=st + i-1)
  xtest <- window(res.fit1, start=st + (i-1) + 1/52, end=st + i)
  
  #models:  ARIMA(0,1,0)x(1,0,0)[52]. and ARIMA(3,0,2)(0,0,0)[52]
  
  fit <- Arima(xtrain, order=c(0,1,0), seasonal=list(order=c(1,0,0), period=52),
                include.drift=TRUE, method="ML")
  fcast <- forecast(fit, h=52)
  
  fit2 <- Arima(xtrain, order=c(3,0,2), seasonal=list(order=c(0,0,0), period=52),
                include.drift=TRUE, method="ML")
  fcast2 <- forecast(fit2, h=52)
  
  

  rmse1[i,1:length(xtest)]  <- sqrt((fcast$mean-xtest)^2)
  rmse2[i,1:length(xtest)] <- sqrt((fcast2$mean-xtest)^2)

  
}

plot(1:52, colMeans(rmse1,na.rm=TRUE), type="l", col=2, xlab="horizon", ylab="RMSE")
lines(1:52, colMeans(rmse2,na.rm=TRUE), type="l",col=3)

legend("topleft",legend=c("fit1: auto.arima model","fit2: manual model"),col=2:3,lty=1)

The auto.arima model had lower RMSE compared to the manual model. Therefore, the best model is ARIMA(0,1,0)(1,0,0)[52].

Code
#canada_value_ts
n=length(res.fit) #205  204
k=48

 #n-k=156; 156/12=13;

res <- res.fit[1:204] 
rmse1 <- matrix(NA, 13,12)
rmse2 <- matrix(NA, 13,12)

st <- tsp(res.fit)[1]+(k-1)/12 

for(i in 1:13)
{
  xtrain <- window(res.fit, end=st + i-1)
  xtest <- window(res.fit, start=st + (i-1) + 1/12, end=st + i)
  
  #models:  ARIMA(1,0,0)x(1,0,0)[12]. and ARIMA(1,0,0)(2,1,1)[12]
  
  fit <- Arima(xtrain, order=c(1,0,0), seasonal=list(order=c(1,0,0), period=12),
                include.drift=TRUE, method="ML")
  fcast <- forecast(fit, h=12)
  
  fit2 <- Arima(xtrain, order=c(1,0,0), seasonal=list(order=c(2,1,1), period=12),
                include.drift=TRUE, method="ML")
  fcast2 <- forecast(fit2, h=12)
  
  

  rmse1[i,1:length(xtest)]  <- sqrt((fcast$mean-xtest)^2)
  rmse2[i,1:length(xtest)] <- sqrt((fcast2$mean-xtest)^2)

  
}

plot(1:12, colMeans(rmse1,na.rm=TRUE), type="l", col=2, xlab="horizon", ylab="RMSE")
lines(1:12, colMeans(rmse2,na.rm=TRUE), type="l",col=3)

legend("topleft",legend=c("fit1: auto.arima model","fit2: manual model"),col=2:3,lty=1)

The manual model had lower RMSE compared to the auto.arima model. Therefore, the best model is ARIMA(1,0,0)(2,1,1)[12].

Model Fitting

Code
#ups_covid_ts
fit <- Arima(ups_covid_ts[,'UPS'], order=c(0,1,0), seasonal=c(1,0,0), xreg=ups_covid_ts[,'COVID'])
summary(fit)
Series: ups_covid_ts[, "UPS"] 
Regression with ARIMA(0,1,0)(1,0,0)[52] errors 

Coefficients:
         sar1     xreg
      -0.0795  -0.0008
s.e.   0.0930   0.0068

sigma^2 = 0.001926:  log likelihood = 289.38
AIC=-572.75   AICc=-572.61   BIC=-563.36

Training set error measures:
                      ME       RMSE        MAE        MPE      MAPE      MASE
Training set 0.003119728 0.04349475 0.03097627 0.06025101 0.6193865 0.1143339
                   ACF1
Training set 0.01248485

GARCH model on SARIMAX residuals

Code
ggAcf((residuals(fit))^2)

Code
ggPacf((residuals(fit))^2)

Code
model <- list() ## set counter
cc <- 1
for (p in 1:6) {
  for (q in 1:6) {
  
model[[cc]] <- garch(residuals(fit),order=c(q,p),trace=F)
cc <- cc + 1
}
} 

## get AIC values for model evaluation
GARCH_AIC <- sapply(model, AIC) ## model with lowest AIC is the best
which(GARCH_AIC == min(GARCH_AIC))
[1] 1
Code
model[[which(GARCH_AIC == min(GARCH_AIC))]]

Call:
garch(x = residuals(fit), order = c(q, p), trace = F)

Coefficient(s):
       a0         a1         b1  
1.805e-03  2.170e-16  4.993e-02  
Code
summary(garchFit(~garch(1,1), residuals(fit),trace = F))

Title:
 GARCH Modelling 

Call:
 garchFit(formula = ~garch(1, 1), data = residuals(fit), trace = F) 

Mean and Variance Equation:
 data ~ garch(1, 1)
<environment: 0x131a9c630>
 [data = residuals(fit)]

Conditional Distribution:
 norm 

Coefficient(s):
        mu       omega      alpha1       beta1  
0.00311950  0.00079769  0.00000001  0.57805135  

Std. Errors:
 based on Hessian 

Error Analysis:
        Estimate  Std. Error  t value Pr(>|t|)    
mu     3.120e-03   3.327e-03    0.938    0.348    
omega  7.977e-04   1.550e-04    5.146 2.67e-07 ***
alpha1 1.000e-08         NaN      NaN      NaN    
beta1  5.781e-01         NaN      NaN      NaN    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log Likelihood:
 292.1893    normalized:  1.71876 

Description:
 Fri Apr 26 00:06:29 2024 by user:  


Standardised Residuals Tests:
                                 Statistic      p-Value
 Jarque-Bera Test   R    Chi^2  79.2303887 0.000000e+00
 Shapiro-Wilk Test  R    W       0.9470052 5.500240e-06
 Ljung-Box Test     R    Q(10)   4.6081684 9.157711e-01
 Ljung-Box Test     R    Q(15)  10.5752144 7.820846e-01
 Ljung-Box Test     R    Q(20)  18.0518484 5.839926e-01
 Ljung-Box Test     R^2  Q(10)   3.7895993 9.563397e-01
 Ljung-Box Test     R^2  Q(15)   6.1393552 9.773106e-01
 Ljung-Box Test     R^2  Q(20)   7.6212699 9.940916e-01
 LM Arch Test       R    TR^2    5.7006263 9.304150e-01

Information Criterion Statistics:
      AIC       BIC       SIC      HQIC 
-3.390462 -3.316678 -3.391536 -3.360521 

Alpha and beta are not significant, so SARIMAX model is enough.

Equation: \(y_t = -0.0008x_t + \Phi(B^{52})(1-B)n_t\)

\(\Phi(B^s) = 1+0.0795B^{52}\)

Code
#canada_value_ts
xreg <- cbind(fuel = canada_value_ts[, "Fuel"],
              tsi = canada_value_ts[, "TSI"])
fit0 <- Arima(canada_value_ts[,'Canada_value'], order=c(1,0,0), seasonal=c(2,1,1), xreg=xreg)
summary(fit0)
Series: canada_value_ts[, "Canada_value"] 
Regression with ARIMA(1,0,0)(2,1,1)[12] errors 

Coefficients:
         ar1    sar1     sar2     sma1    fuel     tsi
      0.7016  0.1968  -0.2513  -0.8508  0.0601  0.3302
s.e.  0.0602  0.0890   0.0850   0.0947  0.0057  0.0557

sigma^2 = 4.181:  log likelihood = -417.99
AIC=849.97   AICc=850.58   BIC=872.81

Training set error measures:
                     ME     RMSE      MAE        MPE     MAPE      MASE
Training set -0.2265331 1.952955 1.423738 -0.6403124 3.005949 0.2500821
                    ACF1
Training set -0.03815489

Equation: \(y_t = 0.0601x_t{fuel} + 0.3302x_t{tsi} + \phi(B) \Phi(B^{12}) (1-B^{12})n_t + \Theta (B^{12}) z_t\)

\(\phi(B) = 1 - 0.7016B\)

\(\Phi(B^s) = 1 - 0.1968B^{12} + 0.2513B^{24}\)

\(\Theta (B^s) = 1 + 0.8508B^{12}\)

Forecast

Code
#ups_covid_ts
covid_fit <- auto.arima(ups_covid_ts[,'COVID'])
fcov <- forecast(covid_fit, 52)

fcast <- forecast(fit, xreg=fcov$mean,52) 
autoplot(fcast, main="Forecast of UPS Stock Price for the next year") + xlab("Year") +
  ylab("UPS Stock Price")

The forecast of the UPS stock price appears to be relatively flat, indicating that the model may not adequately capture the fluctuations observed in the original data.

Code
#canada_value_ts
fuel_fit <- auto.arima(canada_value_ts[,'Fuel'])
ffuel <- forecast(fuel_fit, 36)
tsi_fit <- auto.arima(canada_value_ts[,'TSI'])
ftsi <- forecast(tsi_fit, 36)

fxreg <- cbind(fuel = ffuel$mean, tsi = ftsi$mean) 
fcast <- forecast(fit0, xreg=fxreg,36) 
autoplot(fcast, main="Forecast of Canada Freight Value for the next three years") + xlab("Year") +
  ylab("Canada Freight Value")

The trend of the Canada freight value is observed to be decreasing, and the model appears to capture the seasonal fluctuations present in the original data.

Findings

For the UPS stock price and COVID-19 new cases, there appears to be a correlation between the two variables, with the UPS stock price showing an increase alongside the rise in COVID-19 new cases. However, since the COVID-19 pandemic ended in 2023, the UPS stock price is predicted to remain relatively flat with a slight decrease. However, the model may not adequately capture the fluctuations observed in the original data.

Regarding the Canada freight value, both the fuel price and TSI exhibit correlations with the freight value, showcasing seasonal fluctuations. The trend of the Canada freight value over the next three years is expected to decrease. The model appears to capture the seasonal fluctuations present in the original data.

VAR

Select Variables

Code
#revenue_ts
autoplot(revenue_ts, facets=TRUE) +
  xlab("Year") + ylab("") +
  ggtitle("Air Cargo(Domestic), Air Cargo(International), Fuel Price and TSI in USA")

I have selected Air Cargo (Domestic), Air Cargo (International), Fuel Price, and TSI as variables for analysis. My objective is to explore the interdependencies among these variables. Specifically, I aim to understand whether there is a mutual influence between domestic and international air cargo, whether fuel prices affect air cargo volumes, and whether the TSI correlates with the other variables. By examining the relationships among these variables, I hope to gain insights into the dynamics of the air cargo industry and its drivers.

Code
#emp_ts
autoplot(emp_ts, facets=TRUE) +
  xlab("Year") + ylab("") +
  ggtitle("Employment, Air Cargo(Domestic) and Air Cargo(International) in USA")

I have selected Air Cargo (Domestic), Air Cargo (International), and Employment as variables for analysis. My goal is to examine the potential mutual influence among these variables. Specifically, I aim to understand whether changes in domestic and international air cargo volumes impact employment levels, and vice versa. This analysis will help me understand the interconnectedness of these factors and their implications for economic activity and employment trends.

Code
#air_ts
autoplot(air_ts, facets=TRUE) +
  xlab("Year") + ylab("") +
  ggtitle("Air GDP, Revenue and Employment in USA")

My objective is to explore the potential interrelationships among these variables. Specifically, I aim to understand how changes in GDP and Revenue within the air industry may influence employment levels, and vice versa. By investigating the relationships among these key indicators, I seek to gain insights into the dynamics of the air industry’s economic performance and its impact on employment trends.

Select Best VAR(p)

Based on the selection criterion, models with lag orders p=3 and p=10 were considered. Upon fitting VAR(1), VAR(3), and VAR(10) models and comparing the results, it was observed that the VAR(10) model will be instability. Meanwhile, the VAR(3) model showed stronger correlations among the variables, indicating its potential as the best-fitting model.

Code
#revenue_ts

VARselect(revenue_ts, lag.max=10, type="both")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
    10      3      3     10 

$criteria
                1          2           3           4           5           6
AIC(n) -1.7914947 -2.0213042 -2.59011953 -2.63469636 -2.64432332 -2.69497335
HQ(n)  -1.6472405 -1.7808805 -2.25352637 -2.20193372 -2.11539121 -2.06987177
SC(n)  -1.4338410 -1.4252145 -1.75559402 -1.56173499 -1.33292609 -1.14514027
FPE(n)  0.1667186  0.1325113  0.07505571  0.07183165  0.07121818  0.06780348
                 7           8          9          10
AIC(n) -2.71729130 -2.77094641 -2.8679699 -3.08122759
HQ(n)  -1.99602025 -1.95350588 -1.9543599 -2.07144811
SC(n)  -0.92902236 -0.74424161 -0.6028292 -0.57765108
FPE(n)  0.06644436  0.06314439  0.0575039  0.04665972
Code
summary(vars::VAR(revenue_ts, p=1, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Domestic, International, Fuel, TSI 
Deterministic variables: both 
Sample size: 240 
Log Likelihood: -1115.26 
Roots of the characteristic polynomial:
0.9713 0.9493 0.6821 0.6821
Call:
vars::VAR(y = revenue_ts, p = 1, type = "both")


Estimation results for equation Domestic: 
========================================= 
Domestic = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       0.6579380  0.0645703  10.189  < 2e-16 ***
International.l1  0.0692311  0.0492711   1.405  0.16132    
Fuel.l1          -0.0002871  0.0001188  -2.416  0.01648 *  
TSI.l1            0.0047247  0.0015653   3.018  0.00282 ** 
const            -0.1637126  0.1388310  -1.179  0.23951    
trend            -0.0004562  0.0002529  -1.804  0.07259 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1109 on 234 degrees of freedom
Multiple R-Squared: 0.7496, Adjusted R-squared: 0.7442 
F-statistic: 140.1 on 5 and 234 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Domestic.l1      -7.306e-02  9.349e-02  -0.781    0.435    
International.l1  7.191e-01  7.134e-02  10.081   <2e-16 ***
Fuel.l1           6.131e-05  1.721e-04   0.356    0.722    
TSI.l1            5.227e-03  2.266e-03   2.306    0.022 *  
const             1.980e-02  2.010e-01   0.098    0.922    
trend             1.287e-04  3.662e-04   0.351    0.726    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1606 on 234 degrees of freedom
Multiple R-Squared: 0.7587, Adjusted R-squared: 0.7536 
F-statistic: 147.2 on 5 and 234 DF,  p-value: < 2.2e-16 


Estimation results for equation Fuel: 
===================================== 
Fuel = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + const + trend 

                 Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       4.66268   11.59496   0.402    0.688    
International.l1 14.15808    8.84766   1.600    0.111    
Fuel.l1           0.95592    0.02134  44.793   <2e-16 ***
TSI.l1           -0.41339    0.28109  -1.471    0.143    
const            21.56674   24.93003   0.865    0.388    
trend             0.02366    0.04542   0.521    0.603    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 19.92 on 234 degrees of freedom
Multiple R-Squared: 0.936,  Adjusted R-squared: 0.9346 
F-statistic: 684.4 on 5 and 234 DF,  p-value: < 2.2e-16 


Estimation results for equation TSI: 
==================================== 
TSI = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       0.9331833  0.8919918   1.046   0.2966    
International.l1 -0.0301676  0.6806438  -0.044   0.9647    
Fuel.l1          -0.0005522  0.0016417  -0.336   0.7369    
TSI.l1            0.9428935  0.0216238  43.604   <2e-16 ***
const             4.9601897  1.9178480   2.586   0.0103 *  
trend             0.0080823  0.0034942   2.313   0.0216 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 1.533 on 234 degrees of freedom
Multiple R-Squared: 0.9845, Adjusted R-squared: 0.9841 
F-statistic:  2963 on 5 and 234 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
              Domestic International     Fuel     TSI
Domestic       0.01231       0.01232  -0.1394 0.02611
International  0.01232       0.02580   0.2995 0.05176
Fuel          -0.13935       0.29948 396.9225 2.41270
TSI            0.02611       0.05176   2.4127 2.34903

Correlation matrix of residuals:
              Domestic International     Fuel     TSI
Domestic       1.00000       0.69107 -0.06304 0.15356
International  0.69107       1.00000  0.09358 0.21023
Fuel          -0.06304       0.09358  1.00000 0.07901
TSI            0.15356       0.21023  0.07901 1.00000
Code
summary(vars::VAR(revenue_ts, p=3, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Domestic, International, Fuel, TSI 
Deterministic variables: both 
Sample size: 238 
Log Likelihood: -982.485 
Roots of the characteristic polynomial:
0.9811 0.9539 0.8896 0.8896 0.7354 0.7354 0.6699 0.6699 0.3692 0.3692 0.3253 0.1988
Call:
vars::VAR(y = revenue_ts, p = 3, type = "both")


Estimation results for equation Domestic: 
========================================= 
Domestic = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Domestic.l1      -8.302e-03  7.359e-02  -0.113   0.9103    
International.l1  2.927e-01  4.848e-02   6.038 6.41e-09 ***
Fuel.l1          -2.777e-04  2.864e-04  -0.970   0.3333    
TSI.l1           -8.575e-03  3.657e-03  -2.345   0.0199 *  
Domestic.l2       1.100e-01  7.515e-02   1.463   0.1448    
International.l2  8.956e-02  5.503e-02   1.627   0.1051    
Fuel.l2           7.788e-04  4.362e-04   1.785   0.0755 .  
TSI.l2            3.346e-03  4.821e-03   0.694   0.4883    
Domestic.l3       8.158e-01  7.347e-02  11.104  < 2e-16 ***
International.l3 -4.075e-01  4.852e-02  -8.397 5.25e-15 ***
Fuel.l3          -6.434e-04  2.985e-04  -2.156   0.0322 *  
TSI.l3            6.501e-03  3.571e-03   1.821   0.0700 .  
const             2.351e-02  1.054e-01   0.223   0.8237    
trend             8.982e-05  1.962e-04   0.458   0.6475    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.08101 on 224 degrees of freedom
Multiple R-Squared: 0.8719, Adjusted R-squared: 0.8645 
F-statistic: 117.3 on 13 and 224 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Domestic.l1      -0.6878297  0.1159465  -5.932 1.12e-08 ***
International.l1  0.8730090  0.0763763  11.430  < 2e-16 ***
Fuel.l1           0.0011297  0.0004513   2.503  0.01301 *  
TSI.l1           -0.0073074  0.0057620  -1.268  0.20604    
Domestic.l2      -0.3393576  0.1183990  -2.866  0.00455 ** 
International.l2  0.1950756  0.0867038   2.250  0.02543 *  
Fuel.l2          -0.0001744  0.0006873  -0.254  0.79994    
TSI.l2            0.0062441  0.0075948   0.822  0.41186    
Domestic.l3       1.1293478  0.1157595   9.756  < 2e-16 ***
International.l3 -0.3427729  0.0764511  -4.484 1.17e-05 ***
Fuel.l3          -0.0009782  0.0004703  -2.080  0.03866 *  
TSI.l3            0.0019086  0.0056258   0.339  0.73474    
const             0.2642114  0.1660924   1.591  0.11308    
trend             0.0006541  0.0003091   2.116  0.03542 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1276 on 224 degrees of freedom
Multiple R-Squared: 0.8506, Adjusted R-squared: 0.8419 
F-statistic:  98.1 on 13 and 224 DF,  p-value: < 2.2e-16 


Estimation results for equation Fuel: 
===================================== 
Fuel = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + const + trend 

                 Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       5.28991   17.64864   0.300   0.7647    
International.l1  6.24824   11.62551   0.537   0.5915    
Fuel.l1           1.13224    0.06869  16.484   <2e-16 ***
TSI.l1            0.59519    0.87705   0.679   0.4981    
Domestic.l2      -6.74732   18.02194  -0.374   0.7085    
International.l2 -3.51815   13.19750  -0.267   0.7900    
Fuel.l2          -0.14966    0.10461  -1.431   0.1539    
TSI.l2           -0.19890    1.15603  -0.172   0.8635    
Domestic.l3      -1.02678   17.62018  -0.058   0.9536    
International.l3 20.93866   11.63689   1.799   0.0733 .  
Fuel.l3          -0.05278    0.07158  -0.737   0.4617    
TSI.l3           -0.99149    0.85633  -1.158   0.2482    
const            36.25988   25.28153   1.434   0.1529    
trend             0.03442    0.04704   0.732   0.4651    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 19.43 on 224 degrees of freedom
Multiple R-Squared: 0.9408, Adjusted R-squared: 0.9374 
F-statistic: 273.8 on 13 and 224 DF,  p-value: < 2.2e-16 


Estimation results for equation TSI: 
==================================== 
TSI = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + const + trend 

                  Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       0.372175   1.338447   0.278 0.781219    
International.l1 -0.666864   0.881662  -0.756 0.450221    
Fuel.l1           0.019510   0.005209   3.745 0.000229 ***
TSI.l1            0.845869   0.066514  12.717  < 2e-16 ***
Domestic.l2      -0.630593   1.366757  -0.461 0.644975    
International.l2  1.397715   1.000879   1.396 0.163950    
Fuel.l2          -0.012709   0.007934  -1.602 0.110581    
TSI.l2           -0.043985   0.087671  -0.502 0.616365    
Domestic.l3       0.705271   1.336288   0.528 0.598172    
International.l3 -0.729092   0.882525  -0.826 0.409602    
Fuel.l3          -0.008716   0.005429  -1.606 0.109777    
TSI.l3            0.150518   0.064943   2.318 0.021367 *  
const             4.735929   1.917313   2.470 0.014255 *  
trend             0.007462   0.003568   2.092 0.037608 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 1.473 on 224 degrees of freedom
Multiple R-Squared: 0.9861, Adjusted R-squared: 0.9852 
F-statistic:  1219 on 13 and 224 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
               Domestic International     Fuel     TSI
Domestic       0.006562      0.005779  -0.1260 0.01196
International  0.005779      0.016288   0.1116 0.02800
Fuel          -0.126009      0.111580 377.3863 1.23452
TSI            0.011962      0.028003   1.2345 2.17053

Correlation matrix of residuals:
              Domestic International     Fuel     TSI
Domestic       1.00000        0.5589 -0.08007 0.10023
International  0.55895        1.0000  0.04500 0.14893
Fuel          -0.08007        0.0450  1.00000 0.04313
TSI            0.10023        0.1489  0.04313 1.00000
Code
summary(vars::VAR(revenue_ts, p=10, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Domestic, International, Fuel, TSI 
Deterministic variables: both 
Sample size: 231 
Log Likelihood: -787.217 
Roots of the characteristic polynomial:
0.9895 0.9895 0.9664 0.9664 0.9573 0.9573 0.9568 0.9568 0.9562 0.9562 0.946 0.946 0.9254 0.9074 0.9074 0.8791 0.8791 0.8471 0.8471 0.8458 0.8458 0.8411 0.8411 0.8284 0.8284 0.827 0.827 0.8105 0.7963 0.7963 0.7884 0.7884 0.7835 0.7835 0.7624 0.7624 0.518 0.4269 0.4269 0.0279
Call:
vars::VAR(y = revenue_ts, p = 10, type = "both")


Estimation results for equation Domestic: 
========================================= 
Domestic = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + Domestic.l4 + International.l4 + Fuel.l4 + TSI.l4 + Domestic.l5 + International.l5 + Fuel.l5 + TSI.l5 + Domestic.l6 + International.l6 + Fuel.l6 + TSI.l6 + Domestic.l7 + International.l7 + Fuel.l7 + TSI.l7 + Domestic.l8 + International.l8 + Fuel.l8 + TSI.l8 + Domestic.l9 + International.l9 + Fuel.l9 + TSI.l9 + Domestic.l10 + International.l10 + Fuel.l10 + TSI.l10 + const + trend 

                    Estimate Std. Error t value Pr(>|t|)    
Domestic.l1        6.744e-02  8.081e-02   0.835 0.405005    
International.l1   3.114e-01  5.516e-02   5.644 5.99e-08 ***
Fuel.l1            1.953e-04  2.750e-04   0.710 0.478497    
TSI.l1            -1.069e-02  3.623e-03  -2.950 0.003583 ** 
Domestic.l2       -6.736e-03  8.072e-02  -0.083 0.933583    
International.l2  -1.880e-02  6.564e-02  -0.286 0.774861    
Fuel.l2            4.754e-04  4.067e-04   1.169 0.243973    
TSI.l2             6.340e-03  4.571e-03   1.387 0.167082    
Domestic.l3        7.875e-01  8.465e-02   9.302  < 2e-16 ***
International.l3  -2.491e-01  6.332e-02  -3.934 0.000117 ***
Fuel.l3           -6.193e-04  4.054e-04  -1.528 0.128228    
TSI.l3            -2.390e-03  4.612e-03  -0.518 0.604929    
Domestic.l4       -8.394e-02  1.026e-01  -0.818 0.414297    
International.l4  -6.442e-02  6.510e-02  -0.990 0.323684    
Fuel.l4           -1.774e-04  4.274e-04  -0.415 0.678567    
TSI.l4             3.052e-03  4.751e-03   0.642 0.521399    
Domestic.l5        1.703e-01  1.025e-01   1.661 0.098308 .  
International.l5   6.475e-02  6.831e-02   0.948 0.344380    
Fuel.l5            9.421e-04  4.325e-04   2.178 0.030631 *  
TSI.l5            -1.930e-04  4.635e-03  -0.042 0.966830    
Domestic.l6       -2.293e-01  1.013e-01  -2.263 0.024756 *  
International.l6  -9.264e-02  6.831e-02  -1.356 0.176633    
Fuel.l6           -3.856e-04  4.415e-04  -0.873 0.383609    
TSI.l6             3.323e-03  4.612e-03   0.720 0.472121    
Domestic.l7        1.421e-03  1.001e-01   0.014 0.988695    
International.l7   1.802e-01  6.692e-02   2.693 0.007724 ** 
Fuel.l7           -5.811e-04  4.509e-04  -1.289 0.199002    
TSI.l7            -6.490e-03  4.650e-03  -1.396 0.164462    
Domestic.l8       -2.136e-01  9.774e-02  -2.186 0.030076 *  
International.l8  -2.850e-02  6.636e-02  -0.429 0.668078    
Fuel.l8            9.173e-05  4.629e-04   0.198 0.843148    
TSI.l8             6.749e-03  4.600e-03   1.467 0.143966    
Domestic.l9        3.831e-01  9.931e-02   3.857 0.000157 ***
International.l9   4.151e-02  5.939e-02   0.699 0.485529    
Fuel.l9           -1.064e-04  4.646e-04  -0.229 0.819135    
TSI.l9            -2.271e-03  4.504e-03  -0.504 0.614721    
Domestic.l10       2.112e-02  9.629e-02   0.219 0.826606    
International.l10 -2.465e-01  5.523e-02  -4.463 1.39e-05 ***
Fuel.l10           1.319e-04  3.168e-04   0.416 0.677610    
TSI.l10            5.729e-03  3.366e-03   1.702 0.090408 .  
const             -4.820e-02  1.107e-01  -0.435 0.663727    
trend              1.541e-04  1.986e-04   0.776 0.438650    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.06928 on 189 degrees of freedom
Multiple R-Squared: 0.9206, Adjusted R-squared: 0.9033 
F-statistic: 53.42 on 41 and 189 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + Domestic.l4 + International.l4 + Fuel.l4 + TSI.l4 + Domestic.l5 + International.l5 + Fuel.l5 + TSI.l5 + Domestic.l6 + International.l6 + Fuel.l6 + TSI.l6 + Domestic.l7 + International.l7 + Fuel.l7 + TSI.l7 + Domestic.l8 + International.l8 + Fuel.l8 + TSI.l8 + Domestic.l9 + International.l9 + Fuel.l9 + TSI.l9 + Domestic.l10 + International.l10 + Fuel.l10 + TSI.l10 + const + trend 

                    Estimate Std. Error t value Pr(>|t|)    
Domestic.l1       -0.4463585  0.1307935  -3.413 0.000787 ***
International.l1   0.8262960  0.0892897   9.254  < 2e-16 ***
Fuel.l1            0.0008619  0.0004451   1.937 0.054289 .  
TSI.l1            -0.0053062  0.0058646  -0.905 0.366723    
Domestic.l2       -0.4520497  0.1306542  -3.460 0.000668 ***
International.l2   0.0534592  0.1062504   0.503 0.615449    
Fuel.l2            0.0003634  0.0006584   0.552 0.581623    
TSI.l2             0.0109643  0.0073990   1.482 0.140041    
Domestic.l3        1.0277830  0.1370199   7.501  2.4e-12 ***
International.l3  -0.3927900  0.1024987  -3.832 0.000173 ***
Fuel.l3           -0.0014267  0.0006561  -2.174 0.030923 *  
TSI.l3            -0.0071427  0.0074653  -0.957 0.339899    
Domestic.l4       -0.4696822  0.1660713  -2.828 0.005186 ** 
International.l4   0.2957525  0.1053759   2.807 0.005531 ** 
Fuel.l4           -0.0002808  0.0006918  -0.406 0.685319    
TSI.l4             0.0025034  0.0076898   0.326 0.745126    
Domestic.l5        0.2029304  0.1658830   1.223 0.222726    
International.l5   0.1257090  0.1105670   1.137 0.256999    
Fuel.l5            0.0013849  0.0007001   1.978 0.049361 *  
TSI.l5            -0.0054813  0.0075016  -0.731 0.465877    
Domestic.l6       -0.0822714  0.1640131  -0.502 0.616523    
International.l6  -0.1722818  0.1105664  -1.558 0.120865    
Fuel.l6           -0.0006126  0.0007147  -0.857 0.392424    
TSI.l6             0.0023793  0.0074645   0.319 0.750263    
Domestic.l7        0.1502589  0.1620855   0.927 0.355091    
International.l7   0.2123675  0.1083157   1.961 0.051391 .  
Fuel.l7           -0.0004555  0.0007298  -0.624 0.533244    
TSI.l7            -0.0067914  0.0075265  -0.902 0.368028    
Domestic.l8       -0.5530242  0.1582100  -3.496 0.000590 ***
International.l8   0.1136964  0.1074117   1.059 0.291174    
Fuel.l8           -0.0004348  0.0007493  -0.580 0.562454    
TSI.l8             0.0059292  0.0074451   0.796 0.426805    
Domestic.l9        0.4561202  0.1607527   2.837 0.005045 ** 
International.l9  -0.2206458  0.0961375  -2.295 0.022825 *  
Fuel.l9            0.0005235  0.0007520   0.696 0.487160    
TSI.l9             0.0049403  0.0072910   0.678 0.498857    
Domestic.l10       0.2498502  0.1558528   1.603 0.110579    
International.l10 -0.0950983  0.0894046  -1.064 0.288829    
Fuel.l10          -0.0002251  0.0005128  -0.439 0.661190    
TSI.l10           -0.0024540  0.0054484  -0.450 0.652937    
const              0.4305504  0.1791452   2.403 0.017212 *  
trend              0.0008722  0.0003215   2.713 0.007279 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1121 on 189 degrees of freedom
Multiple R-Squared: 0.8935, Adjusted R-squared: 0.8704 
F-statistic: 38.67 on 41 and 189 DF,  p-value: < 2.2e-16 


Estimation results for equation Fuel: 
===================================== 
Fuel = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + Domestic.l4 + International.l4 + Fuel.l4 + TSI.l4 + Domestic.l5 + International.l5 + Fuel.l5 + TSI.l5 + Domestic.l6 + International.l6 + Fuel.l6 + TSI.l6 + Domestic.l7 + International.l7 + Fuel.l7 + TSI.l7 + Domestic.l8 + International.l8 + Fuel.l8 + TSI.l8 + Domestic.l9 + International.l9 + Fuel.l9 + TSI.l9 + Domestic.l10 + International.l10 + Fuel.l10 + TSI.l10 + const + trend 

                    Estimate Std. Error t value Pr(>|t|)    
Domestic.l1         0.033277  21.769033   0.002   0.9988    
International.l1   14.868425  14.861219   1.000   0.3184    
Fuel.l1             1.047580   0.074077  14.142   <2e-16 ***
TSI.l1              1.013645   0.976085   1.038   0.3004    
Domestic.l2         3.493517  21.745836   0.161   0.8725    
International.l2  -14.833976  17.684117  -0.839   0.4026    
Fuel.l2            -0.062276   0.109576  -0.568   0.5705    
TSI.l2              0.771812   1.231469   0.627   0.5316    
Domestic.l3         0.787114  22.805336   0.035   0.9725    
International.l3    5.510652  17.059697   0.323   0.7470    
Fuel.l3            -0.222845   0.109208  -2.041   0.0427 *  
TSI.l3             -1.131478   1.242513  -0.911   0.3636    
Domestic.l4       -24.676248  27.640593  -0.893   0.3731    
International.l4    5.586249  17.538574   0.319   0.7504    
Fuel.l4             0.130721   0.115145   1.135   0.2577    
TSI.l4              1.237665   1.279870   0.967   0.3348    
Domestic.l5         6.793340  27.609258   0.246   0.8059    
International.l5    7.762524  18.402565   0.422   0.6736    
Fuel.l5            -0.140389   0.116519  -1.205   0.2298    
TSI.l5              1.017841   1.248556   0.815   0.4160    
Domestic.l6        39.617148  27.298042   1.451   0.1484    
International.l6  -10.057402  18.402470  -0.547   0.5854    
Fuel.l6             0.058812   0.118947   0.494   0.6216    
TSI.l6             -1.830923   1.242370  -1.474   0.1422    
Domestic.l7         9.822888  26.977219   0.364   0.7162    
International.l7  -18.863814  18.027861  -1.046   0.2967    
Fuel.l7             0.097762   0.121462   0.805   0.4219    
TSI.l7             -0.269893   1.252692  -0.215   0.8296    
Domestic.l8       -62.219554  26.332181  -2.363   0.0191 *  
International.l8   24.782647  17.877400   1.386   0.1673    
Fuel.l8            -0.222809   0.124716  -1.787   0.0756 .  
TSI.l8             -1.547263   1.239150  -1.249   0.2133    
Domestic.l9         3.539384  26.755380   0.132   0.8949    
International.l9   16.804203  16.000948   1.050   0.2950    
Fuel.l9             0.226140   0.125153   1.807   0.0724 .  
TSI.l9             -0.282481   1.213495  -0.233   0.8162    
Domestic.l10       32.487950  25.939859   1.252   0.2120    
International.l10  18.481716  14.880344   1.242   0.2158    
Fuel.l10           -0.007285   0.085352  -0.085   0.9321    
TSI.l10            -0.175500   0.906816  -0.194   0.8467    
const              43.992265  29.816588   1.475   0.1418    
trend               0.029595   0.053502   0.553   0.5808    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 18.66 on 189 degrees of freedom
Multiple R-Squared: 0.9497, Adjusted R-squared: 0.9387 
F-statistic: 86.96 on 41 and 189 DF,  p-value: < 2.2e-16 


Estimation results for equation TSI: 
==================================== 
TSI = Domestic.l1 + International.l1 + Fuel.l1 + TSI.l1 + Domestic.l2 + International.l2 + Fuel.l2 + TSI.l2 + Domestic.l3 + International.l3 + Fuel.l3 + TSI.l3 + Domestic.l4 + International.l4 + Fuel.l4 + TSI.l4 + Domestic.l5 + International.l5 + Fuel.l5 + TSI.l5 + Domestic.l6 + International.l6 + Fuel.l6 + TSI.l6 + Domestic.l7 + International.l7 + Fuel.l7 + TSI.l7 + Domestic.l8 + International.l8 + Fuel.l8 + TSI.l8 + Domestic.l9 + International.l9 + Fuel.l9 + TSI.l9 + Domestic.l10 + International.l10 + Fuel.l10 + TSI.l10 + const + trend 

                    Estimate Std. Error t value Pr(>|t|)    
Domestic.l1        0.4073345  1.6476679   0.247 0.805007    
International.l1   1.1347232  1.1248250   1.009 0.314360    
Fuel.l1            0.0192443  0.0056068   3.432 0.000735 ***
TSI.l1             0.7460737  0.0738785  10.099  < 2e-16 ***
Domestic.l2        0.2397577  1.6459121   0.146 0.884338    
International.l2   1.4537703  1.3384863   1.086 0.278805    
Fuel.l2           -0.0098351  0.0082937  -1.186 0.237169    
TSI.l2            -0.0561804  0.0932082  -0.603 0.547403    
Domestic.l3        4.1652781  1.7261042   2.413 0.016772 *  
International.l3  -2.7830348  1.2912248  -2.155 0.032399 *  
Fuel.l3           -0.0027994  0.0082658  -0.339 0.735233    
TSI.l3             0.2082287  0.0940441   2.214 0.028013 *  
Domestic.l4       -0.3370996  2.0920781  -0.161 0.872162    
International.l4   1.5148130  1.3274703   1.141 0.255261    
Fuel.l4           -0.0066844  0.0087152  -0.767 0.444050    
TSI.l4             0.0059001  0.0968716   0.061 0.951498    
Domestic.l5       -0.4312724  2.0897064  -0.206 0.836717    
International.l5  -0.1437425  1.3928646  -0.103 0.917914    
Fuel.l5            0.0159222  0.0088192   1.805 0.072602 .  
TSI.l5             0.0947401  0.0945015   1.003 0.317372    
Domestic.l6        3.0830586  2.0661509   1.492 0.137321    
International.l6  -3.7565972  1.3928574  -2.697 0.007628 ** 
Fuel.l6           -0.0215666  0.0090029  -2.396 0.017573 *  
TSI.l6             0.1219059  0.0940333   1.296 0.196414    
Domestic.l7       -2.3924532  2.0418682  -1.172 0.242793    
International.l7   1.4672328  1.3645037   1.075 0.283617    
Fuel.l7            0.0023095  0.0091933   0.251 0.801924    
TSI.l7            -0.0349906  0.0948145  -0.369 0.712509    
Domestic.l8       -4.2118813  1.9930462  -2.113 0.035888 *  
International.l8   0.3175975  1.3531156   0.235 0.814683    
Fuel.l8           -0.0012144  0.0094396  -0.129 0.897772    
TSI.l8            -0.1361725  0.0937895  -1.452 0.148189    
Domestic.l9        1.1229519  2.0250776   0.555 0.579877    
International.l9  -0.7172125  1.2110895  -0.592 0.554422    
Fuel.l9            0.0043073  0.0094727   0.455 0.649843    
TSI.l9             0.0378629  0.0918477   0.412 0.680634    
Domestic.l10      -1.9506506  1.9633519  -0.994 0.321722    
International.l10  1.1596651  1.1262726   1.030 0.304491    
Fuel.l10          -0.0003438  0.0064601  -0.053 0.957619    
TSI.l10           -0.0239242  0.0686357  -0.349 0.727802    
const              4.8281201  2.2567762   2.139 0.033687 *  
trend              0.0067474  0.0040495   1.666 0.097323 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 1.413 on 189 degrees of freedom
Multiple R-Squared: 0.9887, Adjusted R-squared: 0.9862 
F-statistic: 401.7 on 41 and 189 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
               Domestic International       Fuel      TSI
Domestic       0.004799      0.004927  -0.027522 0.006024
International  0.004927      0.012574  -0.002224 0.028267
Fuel          -0.027522     -0.002224 348.316937 2.527075
TSI            0.006024      0.028267   2.527075 1.995426

Correlation matrix of residuals:
              Domestic International      Fuel     TSI
Domestic       1.00000      0.634273 -0.021287 0.06155
International  0.63427      1.000000 -0.001063 0.17845
Fuel          -0.02129     -0.001063  1.000000 0.09585
TSI            0.06155      0.178454  0.095855 1.00000

Based on the selection criterion, models with lag orders p=3, p=4 and p=10 were considered. Upon fitting VAR(1), VAR(3), and VAR(4) models and comparing the results, it was observed that the VAR(10) model will be instability. Meanwhile, the VAR(3) model showed stronger correlations among the variables, indicating its potential as the best-fitting model.

Code
#emp_ts


VARselect(emp_ts, lag.max=10, type="both")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
    10      4      3     10 

$criteria
                 1           2           3           4           5           6
AIC(n)    9.110946    8.773224    7.963614    7.900765    7.930169    7.954179
HQ(n)     9.208607    8.929482    8.178469    8.174216    8.262217    8.344823
SC(n)     9.352448    9.159626    8.494917    8.576969    8.751274    8.920185
FPE(n) 9054.112803 6459.709303 2875.307462 2701.022106 2782.990414 2852.620477
                 7           8           9          10
AIC(n)    7.989814    7.939700    7.745777    7.566330
HQ(n)     8.439055    8.447538    8.312211    8.191361
SC(n)     9.100720    9.195508    9.146485    9.111939
FPE(n) 2958.915419 2817.771246 2324.690271 1946.575388
Code
summary(vars::VAR(emp_ts, p=1, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Employment, International, Domestic 
Deterministic variables: both 
Sample size: 216 
Log Likelihood: -1889.312 
Roots of the characteristic polynomial:
0.9571 0.7598 0.675
Call:
vars::VAR(y = emp_ts, p = 1, type = "both")


Estimation results for equation Employment: 
=========================================== 
Employment = Employment.l1 + International.l1 + Domestic.l1 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1     9.573e-01  2.031e-02  47.127   <2e-16 ***
International.l1  5.940e+03  3.065e+03   1.938   0.0539 .  
Domestic.l1      -2.827e+03  3.993e+03  -0.708   0.4797    
const             1.174e+04  9.453e+03   1.242   0.2155    
trend            -3.231e+00  1.004e+01  -0.322   0.7479    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 7005 on 211 degrees of freedom
Multiple R-Squared: 0.923,  Adjusted R-squared: 0.9216 
F-statistic: 632.5 on 4 and 211 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Employment.l1 + International.l1 + Domestic.l1 + const + trend 

                  Estimate Std. Error t value Pr(>|t|)    
Employment.l1    1.304e-07  4.809e-07   0.271  0.78648    
International.l1 7.150e-01  7.255e-02   9.856  < 2e-16 ***
Domestic.l1      3.018e-02  9.452e-02   0.319  0.74984    
const            4.149e-01  2.238e-01   1.854  0.06515 .  
trend            7.727e-04  2.377e-04   3.251  0.00134 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1658 on 211 degrees of freedom
Multiple R-Squared: 0.7144, Adjusted R-squared: 0.709 
F-statistic: 131.9 on 4 and 211 DF,  p-value: < 2.2e-16 


Estimation results for equation Domestic: 
========================================= 
Domestic = Employment.l1 + International.l1 + Domestic.l1 + const + trend 

                  Estimate Std. Error t value Pr(>|t|)    
Employment.l1    3.365e-07  3.295e-07   1.021   0.3084    
International.l1 6.382e-02  4.972e-02   1.284   0.2007    
Domestic.l1      7.196e-01  6.477e-02  11.109   <2e-16 ***
const            1.107e-02  1.534e-01   0.072   0.9425    
trend            3.914e-04  1.629e-04   2.403   0.0171 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1136 on 211 degrees of freedom
Multiple R-Squared: 0.7541, Adjusted R-squared: 0.7495 
F-statistic: 161.8 on 4 and 211 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
              Employment International Domestic
Employment     4.907e+07     137.91323 19.57489
International  1.379e+02       0.02750  0.01333
Domestic       1.957e+01       0.01333  0.01291

Correlation matrix of residuals:
              Employment International Domestic
Employment       1.00000        0.1187  0.02459
International    0.11873        1.0000  0.70741
Domestic         0.02459        0.7074  1.00000
Code
summary(vars::VAR(emp_ts, p=3, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Employment, International, Domestic 
Deterministic variables: both 
Sample size: 214 
Log Likelihood: -1736.085 
Roots of the characteristic polynomial:
0.975 0.8715 0.8715 0.7398 0.7398 0.6188 0.6188 0.5663 0.5663
Call:
vars::VAR(y = emp_ts, p = 3, type = "both")


Estimation results for equation Employment: 
=========================================== 
Employment = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1     1.425e+00  6.803e-02  20.944  < 2e-16 ***
International.l1  1.666e+03  3.721e+03   0.448   0.6548    
Domestic.l1      -1.269e+04  5.778e+03  -2.197   0.0292 *  
Employment.l2    -7.481e-01  1.097e-01  -6.817 1.04e-10 ***
International.l2  8.677e+03  4.274e+03   2.030   0.0437 *  
Domestic.l2       1.520e+03  5.926e+03   0.256   0.7978    
Employment.l3     2.702e-01  6.791e-02   3.978 9.67e-05 ***
International.l3 -3.540e+03  3.660e+03  -0.967   0.3346    
Domestic.l3       9.636e+03  5.614e+03   1.716   0.0876 .  
const             1.385e+04  8.805e+03   1.573   0.1172    
trend            -8.734e+00  9.378e+00  -0.931   0.3527    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 6140 on 203 degrees of freedom
Multiple R-Squared: 0.9423, Adjusted R-squared: 0.9395 
F-statistic: 331.6 on 10 and 203 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1    -4.255e-06  1.407e-06  -3.024  0.00281 ** 
International.l1  9.449e-01  7.694e-02  12.280  < 2e-16 ***
Domestic.l1      -7.559e-01  1.195e-01  -6.326 1.58e-09 ***
Employment.l2     4.640e-06  2.269e-06   2.045  0.04218 *  
International.l2  2.152e-01  8.840e-02   2.435  0.01577 *  
Domestic.l2      -3.569e-01  1.226e-01  -2.912  0.00399 ** 
Employment.l3    -7.567e-07  1.405e-06  -0.539  0.59067    
International.l3 -4.129e-01  7.569e-02  -5.455 1.42e-07 ***
Domestic.l3       1.302e+00  1.161e-01  11.216  < 2e-16 ***
const             4.174e-01  1.821e-01   2.293  0.02290 *  
trend             5.139e-04  1.939e-04   2.650  0.00869 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.127 on 203 degrees of freedom
Multiple R-Squared: 0.8364, Adjusted R-squared: 0.8283 
F-statistic: 103.8 on 10 and 203 DF,  p-value: < 2.2e-16 


Estimation results for equation Domestic: 
========================================= 
Domestic = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1    -3.145e-06  8.588e-07  -3.662 0.000319 ***
International.l1  3.286e-01  4.697e-02   6.996 3.75e-11 ***
Domestic.l1      -4.441e-02  7.294e-02  -0.609 0.543249    
Employment.l2     2.900e-06  1.385e-06   2.093 0.037580 *  
International.l2  1.306e-01  5.396e-02   2.421 0.016346 *  
Domestic.l2       2.576e-02  7.481e-02   0.344 0.730964    
Employment.l3     3.638e-07  8.573e-07   0.424 0.671729    
International.l3 -4.269e-01  4.620e-02  -9.239  < 2e-16 ***
Domestic.l3       9.043e-01  7.087e-02  12.759  < 2e-16 ***
const            -1.861e-02  1.112e-01  -0.167 0.867217    
trend             3.228e-04  1.184e-04   2.727 0.006953 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.07752 on 203 degrees of freedom
Multiple R-Squared: 0.8894, Adjusted R-squared: 0.8839 
F-statistic: 163.2 on 10 and 203 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
              Employment International   Domestic
Employment     3.771e+07     75.964305 -24.526455
International  7.596e+01      0.016127   0.005081
Domestic      -2.453e+01      0.005081   0.006009

Correlation matrix of residuals:
              Employment International Domestic
Employment       1.00000       0.09742 -0.05153
International    0.09742       1.00000  0.51617
Domestic        -0.05153       0.51617  1.00000
Code
summary(vars::VAR(emp_ts, p=4, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: Employment, International, Domestic 
Deterministic variables: both 
Sample size: 213 
Log Likelihood: -1706.444 
Roots of the characteristic polynomial:
0.9763 0.9302 0.9302 0.761 0.761 0.6724 0.6713 0.6713 0.6339 0.6339 0.1615 0.1615
Call:
vars::VAR(y = emp_ts, p = 4, type = "both")


Estimation results for equation Employment: 
=========================================== 
Employment = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + Employment.l4 + International.l4 + Domestic.l4 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1     1.411e+00  7.090e-02  19.895  < 2e-16 ***
International.l1  4.118e+03  4.002e+03   1.029   0.3048    
Domestic.l1      -1.208e+04  6.574e+03  -1.837   0.0677 .  
Employment.l2    -6.834e-01  1.232e-01  -5.548 9.12e-08 ***
International.l2  6.302e+03  4.907e+03   1.284   0.2005    
Domestic.l2       4.196e+03  6.478e+03   0.648   0.5179    
Employment.l3     1.685e-01  1.233e-01   1.367   0.1731    
International.l3 -3.935e+03  4.389e+03  -0.896   0.3711    
Domestic.l3       1.348e+04  6.116e+03   2.204   0.0287 *  
Employment.l4     5.606e-02  7.074e-02   0.793   0.4290    
International.l4  1.482e+03  4.371e+03   0.339   0.7349    
Domestic.l4      -9.306e+03  7.856e+03  -1.185   0.2376    
const             1.159e+04  9.046e+03   1.281   0.2015    
trend            -7.188e+00  9.658e+00  -0.744   0.4576    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 6124 on 199 degrees of freedom
Multiple R-Squared: 0.9433, Adjusted R-squared: 0.9396 
F-statistic: 254.5 on 13 and 199 DF,  p-value: < 2.2e-16 


Estimation results for equation International: 
============================================== 
International = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + Employment.l4 + International.l4 + Domestic.l4 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1    -3.735e-06  1.432e-06  -2.609  0.00976 ** 
International.l1  8.954e-01  8.081e-02  11.080  < 2e-16 ***
Domestic.l1      -6.635e-01  1.327e-01  -4.998 1.26e-06 ***
Employment.l2     2.767e-06  2.487e-06   1.113  0.26721    
International.l2  2.450e-01  9.907e-02   2.473  0.01424 *  
Domestic.l2      -3.782e-01  1.308e-01  -2.892  0.00426 ** 
Employment.l3     2.250e-06  2.489e-06   0.904  0.36710    
International.l3 -5.458e-01  8.862e-02  -6.159 3.99e-09 ***
Domestic.l3       1.313e+00  1.235e-01  10.636  < 2e-16 ***
Employment.l4    -1.918e-06  1.428e-06  -1.343  0.18085    
International.l4  2.212e-01  8.825e-02   2.506  0.01300 *  
Domestic.l4      -1.314e-01  1.586e-01  -0.828  0.40843    
const             4.714e-01  1.827e-01   2.581  0.01058 *  
trend             4.101e-04  1.950e-04   2.103  0.03673 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.1236 on 199 degrees of freedom
Multiple R-Squared: 0.8479, Adjusted R-squared: 0.838 
F-statistic: 85.36 on 13 and 199 DF,  p-value: < 2.2e-16 


Estimation results for equation Domestic: 
========================================= 
Domestic = Employment.l1 + International.l1 + Domestic.l1 + Employment.l2 + International.l2 + Domestic.l2 + Employment.l3 + International.l3 + Domestic.l3 + Employment.l4 + International.l4 + Domestic.l4 + const + trend 

                   Estimate Std. Error t value Pr(>|t|)    
Employment.l1    -3.454e-06  8.953e-07  -3.858 0.000154 ***
International.l1  3.206e-01  5.054e-02   6.343 1.48e-09 ***
Domestic.l1      -1.123e-01  8.301e-02  -1.353 0.177567    
Employment.l2     3.399e-06  1.555e-06   2.186 0.030002 *  
International.l2  1.637e-01  6.196e-02   2.642 0.008898 ** 
Domestic.l2      -1.670e-02  8.179e-02  -0.204 0.838407    
Employment.l3    -6.071e-07  1.556e-06  -0.390 0.696903    
International.l3 -3.653e-01  5.542e-02  -6.592 3.82e-10 ***
Domestic.l3       8.588e-01  7.723e-02  11.120  < 2e-16 ***
Employment.l4     8.367e-07  8.932e-07   0.937 0.350007    
International.l4 -1.026e-01  5.519e-02  -1.859 0.064572 .  
Domestic.l4       1.693e-01  9.920e-02   1.706 0.089486 .  
const            -3.298e-02  1.142e-01  -0.289 0.773120    
trend             3.654e-04  1.220e-04   2.997 0.003078 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 0.07733 on 199 degrees of freedom
Multiple R-Squared: 0.8921, Adjusted R-squared: 0.8851 
F-statistic: 126.6 on 13 and 199 DF,  p-value: < 2.2e-16 



Covariance matrix of residuals:
              Employment International  Domestic
Employment     3.750e+07      74.24409 -20.44733
International  7.424e+01       0.01529   0.00550
Domestic      -2.045e+01       0.00550   0.00598

Correlation matrix of residuals:
              Employment International Domestic
Employment       1.00000       0.09805 -0.04318
International    0.09805       1.00000  0.57527
Domestic        -0.04318       0.57527  1.00000

Based on the selection criterion, models with lag orders p=1, and p=2 were considered. Upon fitting VAR(1), and VAR(2) models and comparing the results, it was observed that the model showed no stronger correlations among the variables, indicating VAR(1) as the best-fitting model.

Code
#air_ts

VARselect(air_ts, lag.max=5, type="both")
$selection
AIC(n)  HQ(n)  SC(n) FPE(n) 
     1      1      1      2 

$criteria
          1    2    3    4    5
AIC(n) -Inf -Inf -Inf -Inf -Inf
HQ(n)  -Inf -Inf -Inf -Inf -Inf
SC(n)  -Inf -Inf -Inf -Inf -Inf
FPE(n)  NaN    0    0    0    0
Code
summary(vars::VAR(air_ts, p=1, type='both'))

VAR Estimation Results:
========================= 
Endogenous variables: GDP, Revenue, Employment 
Deterministic variables: both 
Sample size: 9 
Log Likelihood: -24.465 
Roots of the characteristic polynomial:
0.9637 0.9637 0.1452
Call:
vars::VAR(y = air_ts, p = 1, type = "both")


Estimation results for equation GDP: 
==================================== 
GDP = GDP.l1 + Revenue.l1 + Employment.l1 + const + trend 

              Estimate Std. Error t value Pr(>|t|)
GDP.l1          0.7536     1.1871   0.635    0.560
Revenue.l1     -2.4272     2.1730  -1.117    0.327
Employment.l1 -33.5100   148.2595  -0.226    0.832
const         557.4928   721.3612   0.773    0.483
trend          -3.2764     7.2385  -0.453    0.674


Residual standard error: 40.1 on 4 degrees of freedom
Multiple R-Squared: 0.3371, Adjusted R-squared: -0.3259 
F-statistic: 0.5085 on 4 and 4 DF,  p-value: 0.7357 


Estimation results for equation Revenue: 
======================================== 
Revenue = GDP.l1 + Revenue.l1 + Employment.l1 + const + trend 

                Estimate Std. Error t value Pr(>|t|)  
GDP.l1          0.006165   0.225900   0.027   0.9795  
Revenue.l1      0.276354   0.413499   0.668   0.5405  
Employment.l1  32.531561  28.211972   1.153   0.3131  
const         -70.219995 137.266254  -0.512   0.6359  
trend          -3.888946   1.377392  -2.823   0.0477 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


Residual standard error: 7.63 on 4 degrees of freedom
Multiple R-Squared: 0.8204, Adjusted R-squared: 0.6409 
F-statistic: 4.569 on 4 and 4 DF,  p-value: 0.08515 


Estimation results for equation Employment: 
=========================================== 
Employment = GDP.l1 + Revenue.l1 + Employment.l1 + const + trend 

               Estimate Std. Error t value Pr(>|t|)
GDP.l1         0.005871   0.008825   0.665    0.542
Revenue.l1    -0.023206   0.016153  -1.437    0.224
Employment.l1 -0.133032   1.102096  -0.121    0.910
const          8.495512   5.362285   1.584    0.188
trend          0.013508   0.053808   0.251    0.814


Residual standard error: 0.2981 on 4 degrees of freedom
Multiple R-Squared: 0.4899, Adjusted R-squared: -0.02014 
F-statistic: 0.9605 on 4 and 4 DF,  p-value: 0.5151 



Covariance matrix of residuals:
              GDP Revenue Employment
GDP        1607.9  227.80   11.90347
Revenue     227.8   58.22    1.55021
Employment   11.9    1.55    0.08885

Correlation matrix of residuals:
              GDP Revenue Employment
GDP        1.0000  0.7445     0.9959
Revenue    0.7445  1.0000     0.6816
Employment 0.9959  0.6816     1.0000
Code
#summary(vars::VAR(air_ts, p=2, type='both'))

Cross Validation

Code
#revenue_ts

n=240
k=48 #12*4

#n-k=192; 192/12=16;

rmse1 <- matrix(NA, 192,4) #here is n-k
rmse2 <- matrix(NA, 192,4)
rmse3 <- matrix(NA,16,12) #here n-k / 12 
year<-c()

ts_obj1 <- head(revenue_ts, -1)

st <- tsp(ts_obj1 )[1]+(k-1)/12 


for(i in 1:16) #here  is (n-k ) / 12
{
  
  xtrain <- window(ts_obj1, end=st + i-1)
  xtest <- window(ts_obj1, start=st + (i-1) + 1/12, end=st + i)
  
  ######## first Model ############
  fit <- vars::VAR(xtrain, p=4, type='both')
  fcast <- predict(fit, newdata=xtest, n.ahead = 12)
  fdo<-fcast$fcst$Domestic
  fin<-fcast$fcst$International
  ffuel<-fcast$fcst$Fuel
  ftsi <-fcast$fcst$TSI
  
  ff<-data.frame(fdo[,1],fin[,1],ffuel[,1],ftsi[,1]) #collecting the forecasts for 4 variables
  
  year<-st + (i-1) + 1/12 #starting year
  
  ####### convert to a time series object ########
  
  ff<-ts(ff,start=c(year,1),frequency = 12)
  
  a = 12*i-11 #going from 1,13, 15, (1,1+12,1+12+12, 1+12+12+12, ...)
  b= 12*i #12, 24, 36
  
  ##### collecting errors ######
  rmse1[c(a:b),]  <-sqrt((ff-xtest)^2)
  
  
  ######## Second Model ############
  fit2 <- vars::VAR(xtrain, p=3, type='both')
  fcast2 <- predict(fit2, newdata=xtest, n.ahead = 12)
  
  fdo<-fcast2$fcst$Domestic
  fin<-fcast2$fcst$International
  ffuel<-fcast2$fcst$Fuel
  ftsi <-fcast2$fcst$TSI
  
  ff2<-data.frame(fdo[,1],fin[,1],ffuel[,1],ftsi[,1])
  
  year<-st + (i-1) + 1/12
  
  ff2<-ts(ff2,start=c(year,1),frequency = 12)
  
  a = 12*i-11
  b= 12*i
  rmse2[c(a:b),]  <-sqrt((ff2-xtest)^2)
}

yr = rep(c(2008:2023),each =12) #year
qr = rep(paste0("M",1:12),16) #monthly

rmse1 = data.frame(yr,qr,rmse1)
names(rmse1) =c("Year", "Month", "Domestic", "International", "Fuel", "TSI")
rmse2 = data.frame(yr,qr,rmse2)
names(rmse2) =c("Year", "Month", "Domestic", "International", "Fuel", "TSI")

ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Domestic),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Domestic),color = "red") +
  labs(
    title = "CV RMSE for Domestic",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = International),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = International),color = "red") +
  labs(
    title = "CV RMSE for International",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Fuel),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Fuel),color = "red") +
  labs(
    title = "CV RMSE for Fuel",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = TSI),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = TSI),color = "red") +
  labs(
    title = "CV RMSE for TSI",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Based on cross-validation results, it was found that the VAR(3) model exhibits a lower RMSE compared to the VAR(4) model. This lower RMSE indicates that the VAR(3) model provides better predictive performance and is therefore considered the superior model in this context.

Code
#emp_ts
n=217
k=48 #12*4

#n-k=168; 168/12=14;

rmse1 <- matrix(NA, 168,3) #here is n-k
rmse2 <- matrix(NA, 168,3)
rmse3 <- matrix(NA,14,12) #here n-k / 12 
year<-c()

ts_obj2 <- head(emp_ts, -1)

st <- tsp(ts_obj2 )[1]+(k-1)/12 


for(i in 1:14) #here  is (n-k ) / 12
{
  
  xtrain <- window(ts_obj2, end=st + i-1)
  xtest <- window(ts_obj2, start=st + (i-1) + 1/12, end=st + i)
  
  ######## first Model ############
  fit <- vars::VAR(xtrain, p=3, type='both')
  fcast <- predict(fit, newdata=xtest, n.ahead = 12)
  fin<-fcast$fcst$International
  fdo<-fcast$fcst$Domestic
  femp<-fcast$fcst$Employment
  
  ff<-data.frame(fin[,1],fdo[,1],femp[,1]) #collecting the forecasts for 3 variables
  
  year<-st + (i-1) + 1/12 #starting year
  
  ####### convert to a time series object ########
  
  ff<-ts(ff,start=c(year,1),frequency = 12)
  
  a = 12*i-11 #going from 1,13, 15, (1,1+12,1+12+12, 1+12+12+12, ...)
  b= 12*i #12, 24, 36
  
  ##### collecting errors ######
  rmse1[c(a:b),]  <-sqrt((ff-xtest)^2)
  
  
  ######## Second Model ############
  fit2 <- vars::VAR(xtrain, p=4, type='both')
  fcast2 <- predict(fit2, newdata=xtest, n.ahead = 12)
  
  fin<-fcast2$fcst$International
  fdo<-fcast2$fcst$Domestic
  femp<-fcast2$fcst$Employment
  ff2<-data.frame(fin[,1],fdo[,1],femp[,1])
  
  year<-st + (i-1) + 1/12
  
  ff2<-ts(ff2,start=c(year,1),frequency = 12)
  
  a = 12*i-11
  b= 12*i
  rmse2[c(a:b),]  <-sqrt((ff2-xtest)^2)
}

yr = rep(c(2010:2023),each =12) #year
qr = rep(paste0("M",1:12),14) #monthly

rmse1 = data.frame(yr,qr,rmse1)
names(rmse1) =c("Year", "Month","International","Domestic","Employment")
rmse2 = data.frame(yr,qr,rmse2)
names(rmse2) =c("Year", "Month","International","Domestic","Employment")

ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = International),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = International),color = "red") +
  labs(
    title = "CV RMSE for International",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Domestic),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Domestic),color = "red") +
  labs(
    title = "CV RMSE for Domestic",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Employment),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Employment),color = "red") +
  labs(
    title = "CV RMSE for Employment",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Based on cross-validation results, it was found that the VAR(3) model exhibits a lower RMSE compared to the VAR(4) model. This lower RMSE indicates that the VAR(4) model provides better predictive performance and is therefore considered the superior model in this context.

Code
#air_ts
n=10 #372
k=1 #12*4

#n-k=9; 

rmse1 <- matrix(NA, 9,3) #here is n-k
rmse2 <- matrix(NA, 9,3)
rmse3 <- matrix(NA,9) #here n-k
year<-c()

ts_obj <- air_ts

st <- tsp(ts_obj )[1] 


for(i in 1:9) #here  is (n-k ) 
{
  
  xtrain <- window(ts_obj, end=st + i-1)
  xtest <- window(ts_obj, start=st + (i-1) + 1, end=st + i)
  
  ######## first Model ############
  fit <- vars::VAR(ts_obj, p=1, type='both')
  fcast <- predict(fit, n.ahead = 3)
  fgdp<-fcast$fcst$GDP
  fre<-fcast$fcst$Revenue
  femp<-fcast$fcst$Employment
  
  ff<-data.frame(fgdp[,1],fre[,1],femp[,1]) #collecting the forecasts for 3 variables
  
  year<-st + (i-1) #starting year
  
  ####### convert to a time series object ########
  
  ff<-ts(ff,start=c(year,1),frequency = 1)
  
  ##### collecting errors ######
  rmse1[i-1,]  <-sqrt((ff-xtest)^2)
  
  
  ######## Second Model ############
  fit2 <- vars::VAR(ts_obj, p=2, type='both')
  fcast2 <- predict(fit2, n.ahead = 3)
  
  fgdp<-fcast2$fcst$GDP
  fre<-fcast2$fcst$Revenue
  femp<-fcast2$fcst$Employment
  ff2<-data.frame(fgdp[,1],fre[,1],femp[,1])
  
  year<-st + (i-1)
  
  ff2<-ts(ff2,start=c(year,1),frequency = 1)
  
  rmse2[i-1,]  <-sqrt((ff2-xtest)^2)
}

yr = rep(c(2013:2021),each =1) #year

rmse1 = data.frame(yr,rmse1)
names(rmse1) =c("Year", "GDP","Revenue","Employment")
rmse2 = data.frame(yr,rmse2)
names(rmse2) =c("Year", "GDP","Revenue","Employment")

ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = GDP),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = GDP),color = "red") +
  labs(
    title = "CV RMSE for GDP",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Revenue),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Revenue),color = "red") +
  labs(
    title = "CV RMSE for Revenue",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Code
ggplot() + 
  geom_line(data = rmse1, aes(x = Year, y = Employment),color = "blue") +
  geom_line(data = rmse2, aes(x = Year, y = Employment),color = "red") +
  labs(
    title = "CV RMSE for Employment",
    x = "Date",
    y = "RMSE",
    guides(colour=guide_legend(title="Fit")))

Based on cross-validation results, it was found that the VAR(1) model exhibits a lower RMSE compared to the VAR(2) model. This lower RMSE indicates that the VAR(1) model provides better predictive performance and is therefore considered the superior model in this context.

Forecast

The most optimal model I fitted is the VAR(3) model.

Code
#revenue_ts
var <- vars::VAR(ts_obj1, p=3, type='both')

forecast(var) %>%
  autoplot() + xlab("Year")

The forecast of the variables appears to capture the trend of the original data adequately. While the fuel trend shows a slight decrease, the trends in other variables appear relatively stable or flat. Furthermore, the model seems to struggle in capturing the seasonality present in the data.

The most optimal model I fitted is the VAR(4) model.

Code
#emp_ts
var <- vars::VAR(ts_obj2, p=4, type='both')

forecast(var) %>%
  autoplot() + xlab("Year")

The forecast of the variables appears to capture the trend of the original data adequately. While the employment trend shows a decreasing trend, the trends in other variables exhibit fluctuations, with an initial increase followed by a subsequent drop. However, it seems that the model struggles to capture the seasonality present in the data.

The most optimal model I fitted is the VAR(1) model.

Code
#air_ts
var <- vars::VAR(ts_obj, p=1, type='both')

forecast(var) %>%
  autoplot() + xlab("Year")

The forecast of the variables does not appear to capture the trend of the original data adequately. All the variables exhibit fluctuations, characterized by an initial increase followed by a subsequent drop, which repeats over time. Additionally, the model diagnostics show no significant correlation among the variables, suggesting that the model performance may not be optimal.

Findings

The first model appears to capture the trend of the original data adequately. While the fuel trend shows a slight decrease, the trends in other variables appear relatively stable or flat. Furthermore, the model seems to struggle in capturing the seasonality present in the data.

The second model appears to capture the trend of the original data adequately. While the employment trend shows a decreasing trend, the trends in other variables exhibit fluctuations, with an initial increase followed by a subsequent drop. However, it seems that the model struggles to capture the seasonality present in the data.

The third model does not appear to capture the trend of the original data adequately. All the variables exhibit fluctuations, characterized by an initial increase followed by a subsequent drop, which repeats over time. Additionally, the model diagnostics show no significant correlation among the variables, suggesting that the model performance may not be optimal.

Back to top